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Abstract 

This paper presents a new approach to determining microstructural enrichment func- 
tions to local fields in non-periodic heterogeneous materials with applications in Par- 
tition of Unity and Hybrid Finite Element schemes. It is based on a concept of 
aperiodic tiling by the Wang tiles, designed to produce microstructures morphologi- 
cally similar to original media and enrichment functions that satisfy the underlying 
governing equations. An appealing feature of this approach is that the enrichment 
functions are defined only on a small set of square tiles and extended to larger do- 
mains by an inexpensive stochastic tiling algorithm, thereby circumventing the need 
for simulations of large heterogeneous bodies. Feasibility of the proposed method- 
ology is demonstrated on local constructions of the stress fields in two-dimensional 
particulate media. 

Keywords: Tiling; Microstructure optimization; Enrichment functions; Partition of 
Unity; Trefftz method; FFT-based solvers 



1. Introduction 

A detailed mechanical analysis of materials with complex microstructure that 
accounts for the full resolution of microstructural heterogeneities by using classical 
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numerical approaches as the Finite Element method has been found computation- 
ally prohibitive [1] . To overcome this, one option consists of modeling a macroscopic 
problem with the help of homogenization techniques based on effective material prop- 
erties [2, 3, 4]. However, this may lead to a considerable loss of information about the 
fine scale behavior, thereby resulting in an inaccurate assessment of microstructural 
effects on macroscopic response. 

An alternative strategy free of this phenomenon, yet computationally tractable, 
proceeds from enrichment based techniques, such as Partition of Unity (PU) [5] and 
Hybrid Finite Element (HFE) methods [6]. In these environments the approxima- 
tion properties of standard finite element spaces is enhanced by introducing a-priori 
knowledge of subscale behaviour into the system through a class of enrichment func- 
tions [7, 8, 9, 10, 11, 12, 13]. These are, however, associated with additional degrees of 
freedom and nested computational demands to solve a particular subscale boundary 
value problem. 

The objective of this work is thus to develop a new algorithm that allows for 
extending the local constructions from computationally tractable samples to entire 
macroscopic domains. The algorithm keeps generated enrichment fields continuous 
across subscale sample boundaries and consistent in terms of statistical properties 
of original and reconstructed material morphologies. It is based on a small number 
of Wang tiles [14, 15, 16] and the stochastic tiling concept introduced by Cohen et 
al. [17]. 

In 1961, Hao Wang introduced a tiling concept involving square tiles of differ- 
ent codes on their edges referred to as Wang tiles [14]. The tiles are connected 
together so that the adjacent edges have the same code and permit computationally 
efficient graphic reproduction of morphological patterns in real time [17, 16, 18, 15]. 
Their desirable aesthetic properties are attributed to the aperiodicity of the tilings, 
whereas the low computational effort results from the use of a small number of tiles 
to compress the entire morphological information [19]. 

Here, we exploit and extend these principles to provide a basis for an efficient 
generation of microstructure-based enrichment functions applicable in PU- or HFE- 
based algorithms for heterogeneous media. In order to meet additional criteria arising 
from such constructions, the Simulated annealing-based optimization [20, 21] is used 
to arrive at optimal tile sets. The performance of the method is illustrated on the 
construction of tile-based stress fields in a particulate two phase composite medium 
with linear elastic phases. 

The paper structure is as follows. The concept of stochastic Wang tiling is de- 
scribed in Section 2. A discussion on the optimization procedure based on prescribed 
statistical descriptors and compatibility of mechanical fields on contiguous tile edges 
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is given in Section 3. Section 4 comprises numerical examples demonstrating the 
performance of the proposed approach. This is complemented with discussion of the 
two strategies to incorporate the enrichments functions into finite element schemes 
in 5. Final remarks on the current developments and future plans are assembled 
in Section 6. Finally, in Appendix A, we present a brief overview of the stress 
analysis algorithm. 

The following nomenclature is used in the text. Scalar quantities are denoted by 
plain letters, e.g. a or A, matrices by bold serif font, e.g. a or A, and tensorial quan- 
tities by bold letters, e.g. a or A. In addition, we adopt the Mandel representation 
of symmetric second- and fourth-order tensors a and A, e.g. [22], so that 





an 




" ^1111 


^1122 


y2Ani2 


a = 


022 


, A = 


^2211 


^2222 


^2^2212 




_A/2ai2_ 




_A/2y4i2n 


"\/2^1222 


^1212 



in the two-dimensional setting. 

2. Aperiodic tilings by sets of Wang tiles 

Consider the two dimensional Euclidean space discretized by a regular square 
grid. Each grid cell contains specific morphological patterns that are compatible on 
contiguous boundaries (Fig. 1(b)). If there are no missing cells inside the domain, 
the discretization is called valid tiling and a single cell is referred to as the Wang 
tile (Fig. 1(a)) [14]. The tiles have different codes on their edges, enumerated here 
by small Greek letters, and are not allowed to rotate during the tiling procedure. 
The number of distinct tiles is limited, though arranged in such a fashion that no 
sub-sequence periodically repeats. The set of all distinct tiles is referred to as the tile 
set. Sets that enable uncountably many, always aperiodic, tilings are called aperiodic 
sets [16]. The assumption of strictly aperiodic sets can be relaxed, though still being 
capable to tile the plane aperiodically, e.g., when utilizing the algorithm due to 
Cohen et al. [17] briefly described in the following section. Note that such tilings are 
substantial generalizations of periodic paving algorithms, which use identical tiles 
formed by unit cells. 

2.1. Tile set setup 

Favorable properties of a tile set to control nested repetitive effects proceed from 
the tile and edge code diversity (Fig. 1(a)). The number of edge codes n1 in the z-th 
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Figure 1: (a) A Wang tile with edge codes {a, l3,j,S}, (b) an example of aperiodic valid tiling 
with highlighted connectivity across south-eastern and north-western edges. 

spatial direction of the Cartesian coordinates may be chosen arbitrarily, while the 
number of tiles must satisfy 



where n'^^ = (njn2)^ is the number of tiles in the complete set and = 2, . . . , \/n^ 

stands for the number of admissible pairs of north-western (NW) edge codes. 

When designing a tile set, one chooses a particular number of edge codes n\ and 
. The complete set of n*^^ tiles is created by mutually permuting the codes. In order 
to tile the plane, the south-eastern edge codes must match those assigned to NW 
edges. Fig. 1(b). Thus, the created tiles are sorted according to NW combinations. 
Finally, a desired number of tiles is chosen using Eq. (1), in such a way that n^"^ 
of tiles is taken from each NW group. The emerging set of tiles is referred to as 
Wn^ /n\-n2 in the sequel; moreover, we denote the relative frequency of occurrence 
of the c-th code in the tile set by gc, see Fig. 2. 
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Figure 2: Tile set W8/2-2 consisting of 8 tiles with 2 vertical {0,7} and 2 horizontal {P,S} edge 
codes with equal frequencies of occurrence Qa = Qf; = Q-y = qs = \ [17]. 

2.2. Stochastic tiling algorithm 

Since there are tiles associated to each NW group, index of the new tile to be 
placed is selected randomly from the set {1, . . . , n^^} with the uniform probability. 




(1) 
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Beforehand, one must select an appropriate NW group compatible with the eastern 
code of a previously placed tile and the southern code of the tile just above the one 
to be placed (outer edges a and 5 of shaded areas in Fig. 1(b)). Aperiodicity of 
the resulting tiling is guaranteed when assuming that the random generator never 
returns a periodic sequence of numbers and that each NW group contains at least 
two distinct tiles [17]. 

3. Designing optimal tiles 

To simplify the exposition, we limit our attention to two-phase composite media 
formed by a matrix phase and equi-sized disks of radius p and the microstructure 
representation built on the Wang tile set W8/2-2.^ The location of the disks within 
the tiles has to be optimized to achieve (i) good approximation of the original mi- 
crostructure in terms of a given morphological descriptor (Section 3.1) and (ii) mi- 
crostructures that guarantee the compatibility of mechanical fields on contiguous 
tile edges, ^ see Section 3.2. Such criteria originate from different perspectives. The 
first goal aims at capturing the dominant spatial features of original media and as 
such it is closely related to microstructure reconstruction or compression algorithms 
developed, e.g. in [23, 24, 25, 26, 27, 28, 29]. The latter criterion ensures that the 
tiling-generated fields comply with the governing differential equations and can be 
used as basis functions in the aforementioned numerical schemes, see Section 5. 

3.1. Statistical properties of the microstructure 

The most common class of statistical descriptors embodies a set of n-point prob- 
ability functions, applicable to generic heterogeneous media [30]. In this paper, the 
focus is on the two-point probability functions, which capture primary phenomena 
as the phase volume fractions, characteristic microstructural length(s), long-range 
order information and (an)isotropy. 

Consider a domain P C occupied by a two-phase heterogeneous material 
discretized by a regular lattice with Nj' x N'^ pixels, indexed by k e /C with 

}C = ImeZ^ : <mi< ^,i = 1,2 

I 2 2 



^The set W8/2-2 has been chosen since it is the simplest one that allows for aperiodic patterns 
in the stochastic sense [17]. Note that most of the steps of the tile set design can be directly 
generalized to arbitrary media and more complex tile sets. 

^ Henceforth, the term "tiling" stands for "valid tiling" exclusively, thereby excluding invalid 
tilings from the consideration. 
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The distribution of individual phases (disks and matrix) within V is quantified by the 
characteristic function x{^)i which equals 1 when k is occupied by the disk phase and 
is otherwise, cf. Fig. 3(a). Assuming a periodic^ ergodic medium, the two-point 
probability function 5*2 : /C — )■ [0, 1] is then defined as [30] 

'^2('<) = ^E^(^)^(Lk + mJJ, (3) 

me/C 

where N'^ = NfN^ and [•J^: denotes the /C-periodic extension. Noticing that (3) 
has the structure of circular correlation, the two-point probability function can be 
efficiently evaluated using Fast Fourier Transform techniques, see e.g. [31]. 

According to its definition, 5*2 (k) quantifies the probability that two arbitrary 
points separated by k will both be located at the disk phase when randomly selected 
from /C. Denoting by (p the disk volume fraction, < < 1, the two-point probability 
function satisfies 5*2(0) = 0. Moreover, S'2(k) ~ 0^ for ||k|| ^ p indicates that the 
medium does not exhibit artificial long-range order effects, cf. Fig. 3(b). 




Figure 3: (a) An example of two-phase medium formed by equilibrium distribution of 1,300 equi- 
sized disks of volume fraction 26.8% and (b) the two-point probability function 6*2; the sample is 
discretized with 1, 000 x 1, 000 pixels and each disk has the radius of 8 pixels. 

The following procedure is adopted to determine the two-point probability func- 
tion for the tile-based microstructure. First, the set W8/2-2 is used to assemble a 



^Note that periodicity is considered here for the sake of computational efficiency. The tiling- 
generated microstructure is always aperiodic. 
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small tiling Vs C consisting of 4 x 4 tiles, Fig. 4(a). Here, each tile appears with 
the same frequency in order to suppress artificial fluctuations in volume fractions. 
The domain Vs is discretized by an N'^'^ x A*";^-^ regular grid with the same lattice 
spacing as in the original microstructure, so that Nf'^ < Nf. 

Given a parameter matrix p G P quantifying positions of individual disks, as dis- 
cussed later on, the tile-based morphology is quantified by the two-point probability 
function 5*2 : P x /C — )■ [0, 1], and its proximity to the target microstructure is given 

by 

1 



/5(P) 



52(p,k; 



(4) 



where N'^^ and KP^ are defined analogously as for the original medium T). 
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(b) Tiling V, 



Figure 4: Valid tilings used in optimization with respect to (a) two-point probability functions and 
(b) stress fields; highligthed vertical edges in (b) correspond to edge set contaning 50 equivalent 
edges of code 5. 



3.2. Compatibility of mechanical fields 

The second, yet more complex, goal is to find the tile set morphology that ensures 
the compatibility of physical fields obtained by the tiling algorithm. In particular, 
we choose the stress cr as the quantity of interest, since it is typically used to pre- 
dict the onset of inelastic behavior. Again, for the computational convenience, the 
optimization is performed for a relatively small representative tiling Dg., Fig. 4(b), 
discretized into a N^'' x N'^" points, periodic on external boundaries. 

Assigning the elastic properties to the matrix and the disks, the stress analysis 
proceeds as follows. First, the distribution of individual phases within the tile set 
T>„ is determined for a parameter vector p. Next, using the algorithm outlined in 
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Appendix A, the tiling is subjected to an average strain field 

E= [En E22 72^12]"" 
that results in the distribution of local stress fields 



(5) 



p,k)= cr;?(p,k) 4^(p,k) y2al^^(p,k) 



k G IC^^. 



(6) 



Here, KP'^ is defined analogously to (2), and the superscript w^*^ denotes the stress 
determined for the i-th loadcase, i = 1,2, 3, obtained by successively setting each 
component of the average strain field E to 1, while the remaining components equal 
to 0. 

Similarly to the original Wang idea, the compatibility of the resulting stress fields 
will be enforced via edges. In particular, for each edge code c = 1,2,. . . ,n'^ with 
= n'i + n\, we identify a set S^, formed by contiguous edges '^c,j^j = 1) • • • > '^c 

of 

identical code c and length i, see Fig. 4(b). The stress incompatibility for the i-th 
loadcase and edge code c is quantified by 



I^fl(p) 



1 ^ E E 

- ^ max {t«(p, Z,,(.))} - min {t«(p, Z,,(.))} 

s=l 



(7) 



where T^*-* G denotes the traction vector acting at determined from cr(*\ 
Zcj (s) G KP'^ gives the coordinate of the s-th pixel at the j-th edge with code c, max 
and min operations are understood component-wise and = Yli Collecting 
the contributions from all loadcases and codes, the objective function becomes 



(8) 



1=1 c=l 



3.3. Micro structure parametrization 

In general, the microstructure representation is based on a given Wang tile set 
consisting of ra* tiles of the edge length £ G N, in which we distribute n'^ disks of radius 



''Notice that 2?a- must contain all admissible combinations of tile pairs from the set Wn'/n^-rij 
sharing all edge codes. For example, for the W8/2-2 tile set, Fig. 2, there are 16 combinations of 
basic tiles sharing the code : {2 - 1, 2 - 2, 2 - 7, 2 - 8, 4 - 1, 4 - 2, 4 - 7, 4 - 8, 5 - 1, 5 - 2, 5 - 
7, 5 - 8, 7 - 1, 7 - 2, 7 - 7, 7 - 8} and therefore = 50 > 16, cf. Fig. 4(b). 
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p. The (i-th disk is represented by a triplet {t^'^\ xf\ x^2'^} , where t''^' G {1, . . . j^*} 
denotes the tile index and x^^^ G {!,...,£} specifies the position of the d-ih disk 
within the tile at the j-th direction. The aforementioned parameter vector p is 
obtained collection of these data: 



f[d] A<1\ Ad 



(9) 

d=l 



This also implicitly defines the parameter space V. 

In an admissible configuration, the disks do not overlap each other or corners of 
tiles being associated with. The first constraint refiects the given feature of the orig- 
inal microstructure, Fig. 3(a), whereas the latter one arises as an artifact intrinsic to 
the edge-based tiling algorithm, e.g., [17]. In addition, to maintain the morpholog- 
ical compatibility, any disk intersecting the edge of a given code is also associated 
to tiles sharing the same code, cf. Fig. 5. To emphasize this, we encode a particular 
configuration as n'^jnj^jj!!^, where nj^ denotes the number of disks intersecting the 
edge of code c. 

3.4- Optimization procedure 

In fact, the goals represented by objective functions (4) and (8) are confiicting. 
Minimizing only with respect to the two-point probability function results in non- 
compatible stress fields, whereas the latter criterion drives the system to a periodic 
distribution of disks on tile edges. To achieve a compromise solution, we introduce 
a combined objective function in the form 

/(p)=w/5(p) + /.(p), (10) 

where w denotes a weighting factor to be estimated by a trial-and-error procedure. 
The minimization of the objective function (10) is performed by the well-established 
Simulated Annealing method [20, 21], extended by a re-annealing phase to escape 
from local extremes, e.g. [32]. 

Given the number of disks n'^ and the target volume fraction 0, we initiate the 
algorithm by determining the number of edge disks nj^ related to the c-th code and 
the tile edge length i. However, this problem is not trivial due to the multiple 
presence of the code-related disks, recall Fig. 5. Thus, we are able to resolve it only 
via a heuristic procedure outlined next. We estimate the total area of disks within 
the tile set as 



[n'' + Y^ {2n'q, - l) < ^^ (H) 



c=l 
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with A"^ denoting the area of a single disk (in square pixels). This corresponds to 
the tile-based volume fraction _ 

(12) 

as possible. In addition, we impose 



which should be as close to the target value 
the condition 



n 



Era 
c=l 



t v^ra r 



(£-2p)2 - 2pii-p) ' ^^^^ 

matching the local volume fractions of disks intersecting and non-intersecting edges 
of the tile set. Thus, given the number of disks attached to codes {n^}^!^, Eqs. (12) 
and (13) implicitly define tile edge lengths, say ii and £2, which should be equal to 
each other for the correct tile set setup. In our case, we sequentially check all values 
{n^}cLi such that > and X]c=i — select the configuration with the 

minimum difference \£i — £2 1-''' On the basis of these data, we randomly generate 
positions of individual disks and assign them to randomly selected tiles, until an 
admissible configuration p is obtained. 




Figure 5: An example of an admissible 14{l-l-2-l} configuration (with = — — 1 and 
Tip = 2 code-related disks) and its modification by disk displacements; disk 9 leaves its parent tile 
4 and randomly enters tiles 1, 3, 6 or 8; disk 13 leaves its parent tile 7 and randomly enters tiles 1, 
3, 5 or 7. 

A single loop of the optimization algorithm involves a sequential selection of a 
disk d = {1, . . . , n^}, and its movement given by 



X 



[d] 



X 



[d] 



+ £{U-l), j = l,2. 



(14) 



repeated until a new admissible configuration p is encountered. In Eq. (14), U 
denotes a random variable with a uniform distribution in the interval [0, 1]. If a disk, 
during its displacement, leaves its parent tile by crossing the edge with code c, it is 
randomly assigned to a tile sharing the same code. Fig. 5. 



^Note that the values of £ and n'^ are kept constant during the optimization process, whereas 
the values of are allowed to change, since disks can move freely between tile interiors and edges. 
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The acceptance of a new solution p is driven by the Metropolis criterion [20] 



exp (MZ^) > u, (15) 

where T denotes the algorithmic temperature, initially set to T^ax and gradually 
reduced by a multiplicative constant < 1 once the loop over all ra*^ disks is com- 
pleted. The entire algorithm terminates after n^ax objective function evaluations. 
Moreover, we keep it restarting when the current temperature T is less than the 
threshold value Tmin.^ Such a step was found beneficial, as the resulting problem is 
highly multi-modal and discontinuous due to the presence of edge-constrained disks. 



4. Results 

The potential of the tile-based representation is demonstrated for the two-phase 
composite medium with the disk volume fraction = 26.8%, Fig. 3, with a moderate 
contrast in the material parameters^ and under the assumption of the plane strain 
state. Four distinct sets W8/2-2, differing in the tile edge length the number of 
total and edge disks n'^jn^}"'' used in the configuration and in the weight factor 
w balancing the geometrical and mechanical compatibility, have been examined, 
see Fig. 3. In particular, our aim to demonstrate that the proposed optimization 
algorithm works well for the tile morphology design and that the tile sets based on 
the specific tilings Vs and Va can be used to represent generic particulate media. 

In Fig. 7, we present the disk configurations and two-point probability functions 
5*2 obtained for the domain T) being tiled by optimized tile sets. We observe that 
both reconstructed functions 5*2 exhibit local peaks exceeding the value of 0^, which 
reveals the presence of characteristic length scales of order I in the tiled medium. 
Nevertheless, the local extremes are notably smaller than the value of corresponding 
to a periodic construction, e.g. [26], and their number can be substantially reduced 
by increasing the edge length Fig. 7(b) . 

Such conclusions are further supported by Fig. 8 showing cross-sections of the 
two-point probability functions in the ki direction for two different values of the 
weighting parameter w. The results demonstrate that for higher values of to, the 



^AU results reported below correspond to the following set of parameters: T^ax = 10 ^ , Tinin = 
10-6, r„,it = (T,,ax/r„,i„) 1/200 and n^^x = 10^"^. 

^We set the Young modulus E — 1 and the Poisson ratio v — 0.125 for the matrix phase, while 
£■ = 10 and v = 0.125 for disks. All values are expressed in consistent units in what follows. 
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(a) 10{l-l-0-0}, = 23.6% 10{l-0-l-l}, = 28.0% 

(b) 17{1-1-1-1}, (f) = 27.9%, 17{1-1-1-1}, (f) = 27.9% 

(c) 27{2-l-l-l}, (p = 26.7% 27{2-l-l-l}, = 26.7% 

(d) 32{2- 1-2-1}, = 26.6% 32{l-2-2-l}, = 26.6% 

^ = 10^ w = 10^ 

Figure 6: Optimized sets W8/2-2 obtained for weighting factors w equal to 10"* and 10^ and for 
configurations with (a) n'^ = 10, £ = 42 px, (b) n'^ = 17, £ = 52 px, (c) n'^ = 27, ^ = 64 px 
and (d) — 38, £ — 74 px; n'^jn^}"^]^ refers to configuration of n"^ disks in total with disks 
intersecting edge c, £ is the tile length and (f> is the reconstructed volume fraction. 

short-range phenomena are captured to a high accuracy and the magnitude of local 
extremes are consistently reduced with the increasing number of disks, albeit at a 
small rate. By decreasing the importance of S'2-based reconstruction. Fig. 8(b), the 
discrepancy between the original and reconstructed medium substantially increases 
at short distances, leading even to inconsistent value of the volume fraction for 10 
disks, and the local peaks become more pronounced as the stress-based criterion 
drives the system towards periodic media. 

Fig. 9 illustrates the ability of the optimization algorithm to achieve self-compatible 
stress fields by comparing the distribution of the normal stresses an obtained for an 
initial and the optimized configuration of disks. Clearly, the initial configuration 
exhibits significant stress jumps at contiguous edges, which are reduced to almost 
identical values by the proposed optimization algorithm. 

To what extent are the stress fields obtained for the specific configuration Do- 
useful for a general tiling procedure? To address this question, we consider the 
reconstructed stress field ct based on the fields found at eight tiles from top rows 
of the tiling. Fig. 10(d), ^ assembled according to the template V^. The local relative 
error associated with such replacement is quantified by 



It is worth noting that we have obtained analogous results for different selections of tiles 1-8. 
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Figure 7: Reconstructed microstructures and two point probability functions ^2 for tile sets with 
(a) n'^ = 10 disks and £ — 42 px and (b) n'^ = 38 disks and £ = 74 px and identical weighting 
factors w = 10^. 



Outcomes of this comparison are shown in Fig. 10 in the form of tihng-based 
microstructures (a), distribution of the normal stress cxn due to the loading by 
macroscopic shear strain (b), their reconstructed counterparts (c) and the spatial 
distribution of the relative error (d). For the microstructure generated from tiles 
with n'^ = 10 disks, we observe that the reconstructed field displays significant dis- 
tributed errors in tile interiors. Similarly to S'2-based criterion, these deviations 
are significantly reduced and become highly localized when increasing the number 
of disks. This claim is further supported by Fig. 11, plotting the evolution of the 
average error 

as a function of the number of disks. For both values ofw, we observe approximately 
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linear convergence with respect to the number of disks n'^. Since the tile design has 
been based on different objective functions, this gives as yet another confirmation of 
their quality. Note that the significant compression has been achieved by the tiling- 
based representation: the original microstructure contains ~ 1, 300 disks, whereas 
the most detailed tile-based representation builds only on 38 disks and is capable of 
producing much larger microstructures. 



5. Enrichment strategies for finite element schemes 

Consider a domain O G M."^ approximated by n'^ finite elements f^e of charac- 
teristic size h, with the heterogeneity distribution quantified by a tiling O D O 
based on an optimized tiled set. Fig. 12. The tiling also defines the "reconstructed" 
microstructure-induced perturbation displacement u*^^^ : (9 — )■ or stress fields 
—*(«) . Q _v. ]^3^ introduced in Eq. (A. 6), associated with the i-th loadcase. The aim 
of this Section is to discuss two approaches of incorporating the reconstructed func- 
tions in the finite element schemes independently on the underlying element mesh. 

The first strategy is based on the partition of unity method, introduced by Melenk 
and Babuska [5], and generalized in numerous aspects later on, e.g. [33, 34]. In this 
framework, we express the approximate solution in the form, 

u''{x) = J2 Nn{x) [a„ + U*(a;)bJ for xeO, (18) 

n=l 



14 
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10 20 30 

s [px] 



Figure 9: Distribution of the normal stress components an due to loading by macroscopic strain 
Ei2 = 1 at edges of the tiling 2?cr, obtained for n'^ = 38 disks with w = 10^. The gray/black 
patterns correspond to initial/optimized morphologies, respectively. 



where is the number of nodes in the finite element mesh, Nn : O — ?■ M denotes the 
standard finite element basis functions, a„ G and G stand for matrices of 
the standard and extended degrees of freedom associated with the n-th node. The 
matrix of enrichment functions U : (9 — )■ M^^^ is defined as, cf. [13], 



U {x) 



u 



u. 



•(1) 



.(2) 
■(2) 



.(3)- 
■(3) 



[X] 



(19) 



The ansatz (18) is then employed in the standard Galerkin procedure, following the 
exposition of Fish and Yuan [8, 9] . 

On the contrary, the Trefftz approach builds on the hybrid stress formulations 
developed by Teixeira de Freitas [6], see also [35] for an overview. It is based on the 
element stress approximation in the form 



(T^ix) 



Se{x)die + S {x)\3e for X ^ Vlf. 
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pSxm 



stands for the standard basis functions of 



where, analogously to Eq. (18), Se G 
the Trefftz method associated with m regular degrees of freedom ag; bg G M? denotes 
the extended degrees of freedom and the enrichment functions are obtained as 
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Figure 10: Assessment of tiling-based distribution of local fields, (a) microstructure obtained by 
tile Va-, distribution of (b) local stress fields cr\^i , (c) reconstructed stress fields ctJ^'' and of the 
(d) reconstruction-based error /|^^ 

Note that the regular and enriched basis functions are selected such that the stress cr^ 
remains self-equlibrated. This is complemented with an independent approximation 
of displacements at the element boundary 

u\x) = Nl{x)al for x G T^, (22) 

involving only regular (macroscopic) shape functions and standard set of degrees 
of freedom. The remainder of the formulation follows from the weak form of the 
equilibrium and compatibility equations, which can be converted to the element 
boundaries by virtue of the divergence theorem, cf. [6, 36]. Note that the particular 
scheme developed in [13] eliminates the extended degrees of freedom be by coupling 
them to the average strain in the matrix phase with the stiffness matrix L™ in the 
following fashion: 

be=|^(n"'(/ S,{x)dx]a,. (23) 
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Figure 11: Relative stress error as a function of the number of disks for different weights w. 

This step has been found to substantially increase the computational efficiency, while 
maintaining a high accuracy of the resulting fields. 



6. Conclusions 

In this work, we have proposed a novel approach to the construction of aperiodic 
local fields in heterogeneous media with potential applications in hybrid FE environ- 
ments. The method is based on the Wang tiling concept that allows us to represent 
complex patterns using a limited set of basic tiles, complemented by the Simulated 
Annealing-based algorithm to arrive at the optimal tile set morphology. On the basis 
on the results obtained for stress analyses in two-dimensional particulate media we 
conjecture that: 

• the proposed method provides a robust tool for compression of disordered mi- 
crostructures and can serve as an efficient microstructure generation algorithm, 

• it allows for aperiodic extensions of local, possibly periodic, fields to substan- 
tially larger domains while maintaining their compatibility, 

• the tiling-based fields can be used to produce microstructure-based enrichment 
functions for generalized Partition of Unity methods or hybrid finite element 
schemes. 

We are fully aware that our conclusions are somewhat provisional, in the sense 
that these are based on a single set of tiles and the specific class of microstructures. 
Extending the current results to general setting is in the focus of our current work. 
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Appendix A. Computation of mechanical fields 

As explained earlier in Section 3.2, our objective is to determine local fields within 
the tiling V„ C subjected to a given macroscopic strain field E G M^^^ under the 
periodic boundary conditions. These follow from the solution of the elastic unit cell 
problem [22, 37] 

£{x) = l{Vu{x) + Vu{xy) , V (t{x)=0, (7{x) = L{x)£{x), (A.l) 

in which cr : V„ — t- M^^^ denotes the symmetric V^-pehodic second-order stress 
tensor field, u : M.'^ designates the displacement field, L : V„ ]^2x2x2x2 

is the symmetric positive-definite fourth-order tensor of materials stiffness and e : 
Per — M^^^ is the "Dcr-periodic strain field satisfying 



1 



e{x)dx = E. (A.2) 
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It is well-known [37, 22] that the solution to the unit cell is characterized by the 
Lippman-Schwinger equation 

e{x) + [ r{x- y)5L{y)e{y) dy = E, (A.3) 

where SL = L — Lq, Lq G ]R2x2x2x2 auxiliary positive definite stiffness tensor 
of a reference medium and the fourth-order tensor T : T)„ 

_^ M2X2X2X2 ig 

related 

to the Green function of the problem (A.l) with L{x) = Lq. It admits a compact 
closed-form expression in the Fourier space, e.g. [37, Section 5.3], and its action can 
be efficiently evaluated by the FFT algorithm. This observation is at the heart of an 
iterative scheme due to Moulinec and Suquet [38], which can be applied to arbitrary 
digitized media. 

In our case, we adopt an accelerated version of the original algorithm based on 
observations due to Zeman et al. [39]. Since the sample is discretized by a regular 
N'^" X N'^" bitmap, it is convenient to project the integral equation onto the space 
of trigonometric polynomials, e.g. [40]. This yields the linear system in the form 

(I + B)e = E, (A.4) 



where e G M^^^i "^^2 " stores the unknown stress values at individual pixels, E G 
]^3xAfj "^xATj is the corresponding matrix of overall strains and matrix B is expressed 
product of several matrices 



B 



X 




/ 5Lnii 

5L2211 5L2222 v^5Liii2 I . (A.5) 
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Here, F G C 1 ^ 2 and implement the forward and the inverse Fourier trans- 
form and, e.g., 5Lii22 G M^i ''^^2 stores the corresponding component of the stiff- 
ness tensor at individual pixels, see [39] for more details. The system (A.4) is solved 
using standard conjugate gradient algorithm. Upon convergence, the distribution of 
the local stress field cr is determined from the solution e by Eq. (A. 1)3. Note that the 
construction of the enrichment functions in Section 5 is based on the perturbation 
fields of displacements and stresses 

u*{x) = u{x) - / u{y) dy, cr*{x) = (t{x) - / (T{y) dy,{A.6) 
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instead of the total values. Both quantities can be efficiently obtained from the 
known strain and stress fields in the Fourier space, e.g., [41]. 
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